1 Data Import and Cleaning

1.1 Stimuli

# import stimuli information
df.stimuli =
  read_csv(paste0(data_dir,
                  "fun-puzzles-exp1/production2_puzzles-test.csv"))
  
df.stimuli %>%
  head()
## # A tibble: 6 × 15
##   puzzle_id collection_name level_name layout level_width level_height
##       <dbl> <chr>                <dbl> <chr>        <dbl>        <dbl>
## 1         0 Microban                67 ['###…           7            8
## 2         1 Microban                33 ['###…           7            7
## 3         2 Microban                26 [' ##…           6            8
## 4         3 Sokogen-990602…         13 ['###…           8            6
## 5         4 Sokogen-990602…         19 ['###…           6            6
## 6         5 Sokogen-990602…         41 ['###…           9            7
## # ℹ 9 more variables: sokobanonline__likes_rate <dbl>,
## #   sokobanonline__solve_rate <dbl>,
## #   sokobanonline__shortest_solution <dbl>, enjoyment_cat <chr>,
## #   difficulty_cat <chr>, shortestPath_cat <chr>, stimuli_set <chr>,
## #   boxes <chr>, start_position <chr>

1.2 Metadata

df.metadata =
  read_csv(paste0(data_dir,
                  "exp1_puzzle-metadata.csv"))

df.metadata %>%
  head()
## # A tibble: 6 × 20
##   puzzle_id collection_name level_name layout level_width level_height
##       <dbl> <chr>                <dbl> <chr>        <dbl>        <dbl>
## 1         0 Microban                67 ['###…           7            8
## 2         1 Microban                33 ['###…           7            7
## 3         2 Microban                26 [' ##…           6            8
## 4         3 Sokogen-990602…         13 ['###…           8            6
## 5         4 Sokogen-990602…         19 ['###…           6            6
## 6         5 Sokogen-990602…         41 ['###…           9            7
## # ℹ 14 more variables: sokobanonline__likes_rate <dbl>,
## #   sokobanonline__solve_rate <dbl>,
## #   sokobanonline__shortest_solution <dbl>, enjoyment_cat <chr>,
## #   difficulty_cat <chr>, shortestpath_cat <chr>, stimuli_set <chr>,
## #   boxes <chr>, start_position <chr>, shortest_solution <dbl>,
## #   astar_solution <chr>, astar_iters <dbl>,
## #   astar_solution_length <dbl>, astar_solution_string <chr>

1.3 Exclusions

df.exclusions =
  read_csv(paste0(data_dir,
                  "fun-puzzles-exp1/production2_exclusions.csv"))

df.exclusions %>%
  count(reasons)
## # A tibble: 10 × 2
##    reasons                                          n
##    <chr>                                        <int>
##  1 Failed comprehension                             2
##  2 Only 1 practice trials                           1
##  3 Only 2 practice trials                           5
##  4 bad trace log                                   38
##  5 comparison: only 7 posttest trials               2
##  6 comparison: only 7 pretest trials                4
##  7 empty trace log                                 12
##  8 test attempts: only 7 trials with trace data     6
##  9 test ratings: only 6 trials                      1
## 10 test ratings: only 7 trials                     15

1.4 Exit Survey Data

df.survey =
  read_csv(paste0(data_dir,
                  "fun-puzzles-exp1/production2_survey.csv"))

df.survey_processed =
  df.survey %>%
  # mutate scales to show labels
  mutate(judgedDifficulty = factor(judgedDifficulty,
                                   levels = c(1, 2, 3, 4, 5),
                                   labels = c("1\nVery Easy", "2", "3", "4", "5\nVery Hard")),
         participantEffort = factor(participantEffort,
                                    levels = c(1, 2, 3, 4, 5),
                                    labels = c("1\nVery Low", "2", "3", "4", "5\nVery High")),
         sokobanFamiliarity = factor(sokobanFamiliarity,
                                     levels = c(0, 1, 2, 3),
                                     labels = c("Never", "A few times", "Several times", "Very familiar")),
         gamingFrequency = factor(gamingFrequency,
                                  levels = c(0, 1, 2, 3, 4),
                                  labels = c("None", "< 1h", "1-5h", "6-10h", "> 10h"))) %>%
  # filter out exclusions
  filter(!gameID %in% df.exclusions$gameID)

df.survey_processed %>%
  head()
## # A tibble: 6 × 27
##   project     experiment       iteration   stim_id    gameID condition
##   <chr>       <chr>            <chr>       <chr>      <chr>  <chr>    
## 1 fun-puzzles fun-puzzles-exp1 production2 6790956f9… 3940-… difficult
## 2 fun-puzzles fun-puzzles-exp1 production2 6790956f9… 7209-… difficult
## 3 fun-puzzles fun-puzzles-exp1 production2 6790956f9… 4733-… difficult
## 4 fun-puzzles fun-puzzles-exp1 production2 6790956f9… 6110-… difficult
## 5 fun-puzzles fun-puzzles-exp1 production2 6790956f9… 0047-… difficult
## 6 fun-puzzles fun-puzzles-exp1 production2 6790956f9… 2839-… difficult
## # ℹ 21 more variables: startInstructionTS <dbl>,
## #   startPracticeTS <dbl>, startPreTS <dbl>, startMainTS <dbl>,
## #   startPostTS <dbl>, startSurveyTS <dbl>,
## #   startComprehensionTS <dbl>, endExperimentTS <dbl>,
## #   participantExplanations <chr>, judgedDifficulty <fct>,
## #   participantEffort <fct>, sokobanFamiliarity <fct>,
## #   gamingFrequency <fct>, participantYears <dbl>, …

1.5 Gameplay Data

8 test trials:

df.game =
  read_csv(paste0(data_dir,
                  "/fun-puzzles-exp1/production2_testTrials.csv"))

df.game_processed = 
  df.game %>%
  # filter out exclusions
  filter(!gameID %in% df.exclusions$gameID) %>%
  # complement puzzle ids
  mutate(puzzle_id = paste(str_sub(collection_name, 1, 7), level_name)) %>%
  # rename cols
  rename(rating = rate_response,
         ratingRT = rate_rt)

df.game_processed %>%
  head()
## # A tibble: 6 × 32
##   project    experiment iteration stim_id gameID condition study_phase
##   <chr>      <chr>      <chr>     <chr>   <chr>  <chr>     <chr>      
## 1 fun-puzzl… fun-puzzl… producti… 679095… 3940-… difficult test       
## 2 fun-puzzl… fun-puzzl… producti… 679095… 3940-… difficult test       
## 3 fun-puzzl… fun-puzzl… producti… 679095… 3940-… difficult test       
## 4 fun-puzzl… fun-puzzl… producti… 679095… 3940-… difficult test       
## 5 fun-puzzl… fun-puzzl… producti… 679095… 3940-… difficult test       
## 6 fun-puzzl… fun-puzzl… producti… 679095… 3940-… difficult test       
## # ℹ 25 more variables: trialNum <dbl>, stimuli_set <chr>,
## #   author_name <chr>, collection_name <chr>, level_name <dbl>,
## #   width <dbl>, height <dbl>, layout <chr>, boxes <chr>,
## #   start_position.x <dbl>, start_position.y <dbl>,
## #   select_response <dbl>, select_rt <dbl>, rating <dbl>,
## #   ratingRT <dbl>, solveDuration <dbl>, solved <dbl>,
## #   boxesSolved <dbl>, attempt_nInputEvents <dbl>, …
# calculate summary ratings
df.summary_ratings =
  df.game_processed %>%
  # bootstrap ci
  group_by(stimuli_set, puzzle_id, condition) %>%
  tidyboot_mean(rating, na.rm = T) %>%
  ungroup() %>%
  # median, sd, etc.
  left_join(summarizer(df.game_processed, target_cols = c(rating), 
                       puzzle_id, condition) %>%
              dplyr::select(-ends_with('mean')),
            by = c("puzzle_id", "condition")) %>%
  rename(rating_meanBoot = empirical_stat, 
         rating_mean = mean,
         rating_ciLo = ci_lower,
         rating_ciHi = ci_upper)

df.summary_ratings %>%
  head()
## # A tibble: 6 × 10
##   stimuli_set puzzle_id   condition     n rating_meanBoot rating_ciLo
##   <chr>       <chr>       <chr>     <int>           <dbl>       <dbl>
## 1 A           Microba 104 difficult    36            9.08        8.51
## 2 A           Microba 104 enjoyable    32            4.28        3.41
## 3 A           Microba 67  difficult    36            8.67        8   
## 4 A           Microba 67  enjoyable    32            4.03        3.23
## 5 A           Sokogen 13  difficult    36            7.25        6.29
## 6 A           Sokogen 13  enjoyable    32            5.59        4.57
## # ℹ 4 more variables: rating_mean <dbl>, rating_ciHi <dbl>,
## #   rating_median <dbl>, rating_sd <dbl>
# add avg ratings to game data df
df.game_processed =
  df.game_processed %>%
  left_join(df.summary_ratings %>%
              dplyr::select(condition, puzzle_id, rating_mean) %>%
  pivot_wider(names_from = condition, values_from = rating_mean)) %>%
  rename(mean_difficult = difficult,
         mean_enjoyable = enjoyable) %>%
  mutate(puzzle_id_sorted = reorder(puzzle_id, mean_difficult))

16 comparison trials:

df.comparison =
  read_csv(paste0(data_dir,
                  "fun-puzzles-exp1/production2_compareTrials.csv"))

df.comparison_processed =
  df.comparison %>%
  # filter out exclustions
  filter(!gameID %in% df.exclusions$gameID) %>%
  # keep only relevant cols
  dplyr::select(gameID:level_name_1) %>%
  # add phase and puzzle info
  mutate(study_phase = factor(study_phase,
                              levels = c("pretest", "posttest")),
         puzzle_id_0 = paste(str_sub(collection_name_0, 1, 7),
                             level_name_0),
         puzzle_id_1 = paste(str_sub(collection_name_1, 1, 7),
                             level_name_1)) %>%
  # keep only relevant cols
  dplyr::select(-c(collection_name_0:level_name_1)) %>%
  # compute response times
  mutate(rt_done2 = rt_done2 - rt_done1,
         rt_choose = rt - rt_done2)

df.comparison_processed %>%
  head()
## # A tibble: 6 × 13
##   gameID    condition stim_id study_phase trial_type trialNum response
##   <chr>     <chr>     <chr>   <fct>       <chr>         <dbl>    <dbl>
## 1 3940-a6a… difficult 679095… pretest     sokoban-c…        1        0
## 2 3940-a6a… difficult 679095… pretest     sokoban-c…        2        0
## 3 3940-a6a… difficult 679095… pretest     sokoban-c…        3        1
## 4 3940-a6a… difficult 679095… pretest     sokoban-c…        4        0
## 5 3940-a6a… difficult 679095… pretest     sokoban-c…        5        0
## 6 3940-a6a… difficult 679095… pretest     sokoban-c…        6        1
## # ℹ 6 more variables: rt <dbl>, rt_done1 <dbl>, rt_done2 <dbl>,
## #   puzzle_id_0 <chr>, puzzle_id_1 <chr>, rt_choose <dbl>

Create dataset containing all relevant information at once:

df.analysis =
  df.game_processed %>%
  left_join(df.survey_processed %>%
              dplyr::select(gameID,
                     judgedDifficulty,
                     participantEffort,
                     sokobanFamiliarity),
            by = "gameID") %>%
  left_join(df.metadata %>%
              dplyr::select(puzzle_id_num = puzzle_id,
                     layout,
                     sokobanonline__solve_rate,
                     astar_iters, astar_solution_length),
            by = "layout") %>%
  left_join(df.survey_processed %>%
            dplyr::select(gameID,
                   age = participantYears,
                   gender = participantGender),
          by = "gameID") %>%
  mutate(judgedDifficulty = as.numeric(judgedDifficulty),
         sokobanFamiliarity = as.numeric(sokobanFamiliarity),
         participantEffort = as.numeric(participantEffort)) %>% 
  dplyr::select(subject_id = gameID,
         age, gender,
         puzzle_id_sorted, puzzle_id_num,
         trialNum, stimuli_set,
         layout,
         select_response, rating,
         solveDuration, solved, boxesSolved,
         attempt_nInputEvents, attempt_nRestart, attempt_nUndo,
         sokobanonline__solve_rate, astar_iters, astar_solution_length,
         judgedDifficulty, mean_difficult,
         participantEffort, mean_enjoyable,
         sokobanFamiliarity)

df.analysis %>%
  head()
## # A tibble: 6 × 24
##   subject_id        age gender puzzle_id_sorted puzzle_id_num trialNum
##   <chr>           <dbl> <chr>  <fct>                    <dbl>    <dbl>
## 1 3940-a6a959f6-…    23 Female Microba 67                   0        1
## 2 3940-a6a959f6-…    23 Female Sokogen 13                   3        2
## 3 3940-a6a959f6-…    23 Female Sokogen 32                   6        3
## 4 3940-a6a959f6-…    23 Female Sokogen 18                   9        4
## 5 3940-a6a959f6-…    23 Female Microba 104                 12        5
## 6 3940-a6a959f6-…    23 Female Sokogen 38                  15        6
## # ℹ 18 more variables: stimuli_set <chr>, layout <chr>,
## #   select_response <dbl>, rating <dbl>, solveDuration <dbl>,
## #   solved <dbl>, boxesSolved <dbl>, attempt_nInputEvents <dbl>,
## #   attempt_nRestart <dbl>, attempt_nUndo <dbl>,
## #   sokobanonline__solve_rate <dbl>, astar_iters <dbl>,
## #   astar_solution_length <dbl>, judgedDifficulty <dbl>,
## #   mean_difficult <dbl>, participantEffort <dbl>, …

2 Descriptives

Age by condition:

df.analysis %>%
  distinct(subject_id, age, .keep_all = F) %>%
  summarise(n = n(),
            M = mean(age),
            SD = sd(age),
            min = min(age),
            max = max(age))
## # A tibble: 1 × 5
##       n     M    SD   min   max
##   <int> <dbl> <dbl> <dbl> <dbl>
## 1   205  36.7  12.3    18    74

Gender:

df.analysis %>%
  distinct(subject_id, gender, .keep_all = F) %>%
  count(gender)
## # A tibble: 4 × 2
##   gender                   n
##   <chr>                <int>
## 1 Female                  87
## 2 Male                   114
## 3 Non-binary               3
## 4 Prefer not to answer     1

Experience:

df.survey_processed %>%
  count(sokobanFamiliarity)
## # A tibble: 4 × 2
##   sokobanFamiliarity     n
##   <fct>              <int>
## 1 Never                126
## 2 A few times           53
## 3 Several times         18
## 4 Very familiar          8

3 Plots

3.1 % solved per puzzle

df.analysis %>%
  group_by(puzzle_id_sorted) %>%
  summarise(n = n(),
            M = mean(solved, na.rm = T),
            SD = sd(solved, na.rm = T),
            min = min(solved, na.rm = T),
            max = max(solved, na.rm = T)) %>%
    arrange(desc(M))
## # A tibble: 24 × 6
##    puzzle_id_sorted     n     M    SD   min   max
##    <fct>            <int> <dbl> <dbl> <dbl> <dbl>
##  1 Dimitri 37          67 0.806 0.398     0     1
##  2 Sokogen 17          67 0.776 0.420     0     1
##  3 Microba 48          67 0.731 0.447     0     1
##  4 Microba 43          70 0.671 0.473     0     1
##  5 Sokogen 14          70 0.671 0.473     0     1
##  6 Sokogen 6           67 0.627 0.487     0     1
##  7 Microba 41          67 0.567 0.499     0     1
##  8 Sokogen 32          68 0.5   0.504     0     1
##  9 Sokogen 13          68 0.471 0.503     0     1
## 10 Microba 26          70 0.443 0.500     0     1
## # ℹ 14 more rows
df.analysis %>%
  filter(!is.na(solved)) %>%
  mutate(puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                        solved,
                                        .fun = mean, .na_rm = TRUE)) %>%
  group_by(puzzle_id_sorted) %>%
  filter(n() > 0) %>%
  ungroup() %>%
  ggplot(aes(x = factor(puzzle_id_sorted),
             y = (solved)*100,
             fill = stimuli_set)) +
  stat_summary(fun = mean,
               geom = "bar",
               width = 0.6) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "errorbar",
               width = 0.2) +
  scale_y_continuous(limits = c(0, 100),
                     breaks = seq(0, 100, 25)) +
  labs(x = "Puzzle name",
       y = "% solved",
       title = "% solved per puzzle",
       subtitle = "N = 205 subjects, 24 puzzles total (n = 67-70 each)",
       fill = "Puzzle block") +
  coord_flip()

3.2 No of boxes solved per puzzle

df.analysis %>%
  group_by(puzzle_id_sorted) %>%
  summarise(n = n(),
            M = mean(boxesSolved, na.rm = T),
            SD = sd(boxesSolved, na.rm = T),
            min = min(boxesSolved, na.rm = T),
            max = max(boxesSolved, na.rm = T)) %>%
    arrange(desc(M))
## # A tibble: 24 × 6
##    puzzle_id_sorted     n     M    SD   min   max
##    <fct>            <int> <dbl> <dbl> <dbl> <dbl>
##  1 Dimitri 37          67  2.78 0.487     1     3
##  2 Sokogen 17          67  2.64 0.711     1     3
##  3 Microba 43          70  2.57 0.672     1     3
##  4 Microba 48          67  2.52 0.859     0     3
##  5 Sokogen 6           67  2.51 0.726     0     3
##  6 Sokogen 14          70  2.49 0.812     0     3
##  7 Sokogen 32          68  2.41 0.674     0     3
##  8 Sokogen 13          68  2.40 0.626     1     3
##  9 Microba 41          67  2.36 0.811     1     3
## 10 Microba 17          70  2.33 0.557     0     3
## # ℹ 14 more rows
df.analysis %>%
  filter(!is.na(solved)) %>%
  mutate(puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                        solved,
                                        .fun = mean, .na_rm = TRUE)) %>%
  group_by(puzzle_id_sorted) %>%
  filter(n() > 0) %>%
  ungroup() %>%
  ggplot(aes(x = factor(puzzle_id_sorted),
             y = boxesSolved)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "pointrange",
               width = 0.2,
               aes(color = stimuli_set)) +
  geom_point(position = position_jitter(height = .05,
                                      width = 0.1),
           alpha = .05,
           size = 2,
           show.legend = FALSE) +
  scale_y_continuous(limits = c(0, 3),
                     breaks = seq(0, 3, 1)) +
  labs(x = "Puzzle name",
       y = "# of boxes solved",
       title = "Number of boxes solved per puzzle (out of 3)",
       subtitle = "N = 205 subjects, 24 puzzles total (n = 67-70 each)",
       color = "Puzzle block") +
  coord_flip()

3.3 No of actions

df.analysis %>%
  group_by(puzzle_id_sorted) %>%
  summarise(n = n(),
            M = mean(attempt_nInputEvents),
            SD = sd(attempt_nInputEvents),
            min = min(attempt_nInputEvents),
            max = max(attempt_nInputEvents)) %>%
    arrange(desc(M))
## # A tibble: 24 × 6
##    puzzle_id_sorted     n     M    SD   min   max
##    <fct>            <int> <dbl> <dbl> <dbl> <dbl>
##  1 Microba 33          67  265.  129.    43   622
##  2 Sokogen 40          70  236.  139.    22   614
##  3 Microba 104         68  228.  144.    10   608
##  4 Microba 67          68  224.  110.    23   473
##  5 Sokogen 41          70  211.  121.     8   537
##  6 Sokogen 46          70  209.  121.    17   511
##  7 Sokogen 38          68  198.  109.    12   475
##  8 Sokogen 39          67  195.  109.    34   429
##  9 Sokogen 19          67  192.  105.    28   558
## 10 Microba 26          70  190.  126.    23   548
## # ℹ 14 more rows
df.analysis %>%
  filter(!is.na(solved)) %>%
  mutate(puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                        solved,
                                        .fun = mean, .na_rm = TRUE)) %>%
  group_by(puzzle_id_sorted) %>%
  filter(n() > 0) %>%
  ungroup() %>%
  ggplot(aes(x = factor(puzzle_id_sorted),
             y = attempt_nInputEvents)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "pointrange",
               width = 0.2,
               aes(color = stimuli_set)) +
  geom_point(position = position_jitter(height = .05,
                                    width = .1),
         alpha = .05,
         size = 2,
         show.legend = FALSE) +
  scale_y_continuous(limits = c(0, 600),
                   breaks = seq(0, 600, 50)) +
  labs(x = "Puzzle name",
       y = "# of moves attempted",
       title = "Number of moves per puzzle",
       subtitle = "N = 205 subjects, 24 puzzles total (n = 67-70 each)",
       color = "Puzzle block") +
  coord_flip()

3.4 Time until solution

df.analysis %>%
  group_by(puzzle_id_sorted) %>%
  summarise(n = n(),
            M = mean(solveDuration, na.rm = T),
            SD = sd(solveDuration, na.rm = T),
            min = min(solveDuration, na.rm = T),
            max = max(solveDuration, na.rm = T)) %>%
    arrange(desc(M))
## # A tibble: 24 × 6
##    puzzle_id_sorted     n     M    SD   min   max
##    <fct>            <int> <dbl> <dbl> <dbl> <dbl>
##  1 Sokogen 46          70  248.  51.7   160   300
##  2 Microba 104         68  212.  57.5   156   280
##  3 Sokogen 35          68  204.  62.4    82   276
##  4 Sokogen 34          68  197.  58.4   103   299
##  5 Sokogen 40          70  191.  97.8   104   297
##  6 Microba 33          67  191.  82.8    39   295
##  7 Microba 41          67  185.  68.6    63   293
##  8 Sokogen 41          70  184.  64.7    95   297
##  9 Sokogen 39          67  182.  70.1    82   288
## 10 Sokogen 38          68  167.  58.6    81   221
## # ℹ 14 more rows
df.analysis %>%
  filter(!is.na(solved)) %>%
  mutate(puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                 solved,
                                 .fun = mean, .na_rm = TRUE)) %>%
  group_by(puzzle_id_sorted) %>%
  filter(n() > 0) %>%
  ungroup() %>%
  ggplot(aes(x = factor(puzzle_id_sorted),
             y = solveDuration)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "pointrange",
               width = 0.2,
               aes(color = stimuli_set)) +
  geom_point(position = position_jitter(height = .05,
                                        width = .1),
             size = 2,
             alpha = .05,
             show.legend = FALSE) +
  scale_y_continuous(limits = c(0, 300),
                     breaks = seq(0, 300, 60)) +
  labs(x = "Puzzle name",
       y = "Time until solution (in sec)",
       title = "Time until solution per puzzle",
       subtitle = "N = 205 subjects, 24 puzzles total (n = 67-70 each)",
       color = "Puzzle block") +
  coord_flip()

4 Modeling

4.1 Improvement over time

4.1.1 Performance on first half of puzzles vs last half of puzzles

Do people generally become better at Sokoban over time?

I.e., do we see differences in performance metrics (% solved, boxes solved, time until solution, number of steps) when we compare the first four puzzles people solved, vs the last four puzzles people solved?

4.1.1.1 % solved

plot1 =
  df.analysis %>%
  mutate(trialNum = case_when(trialNum > 4 ~ "Last half of puzzles",
                              trialNum <= 4 ~ "First half of puzzles"),
         puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                        solved,
                                        .fun = mean, .na_rm = TRUE)) %>%
  filter(!is.na(trialNum)) %>%
  ggplot(aes(x = puzzle_id_sorted,
             y = (solved)*100,
             fill = trialNum)) +
  stat_summary(fun = mean,
               geom = "bar",
               width = .6,
               position = position_dodge(.8)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "errorbar",
               position = position_dodge(.8)) +
  scale_y_continuous(limits = c(0, 100),
                     breaks = seq(0, 100, 25)) +
  labs(x = "Puzzle name",
       y = "% solved",
       title = "% solved over time",
       fill = "Stage") +
  coord_flip() +
  theme(legend.position = "top")

plot1

Do we see a progression in % solved when we break it down by puzzle order?

df.analysis %>%
  group_by(trialNum, puzzle_id_sorted) %>%
  summarise(mean = mean(solved)) %>%
  ggplot(aes(x = trialNum,
             y = (mean*100),
             color = puzzle_id_sorted)) +
  geom_point(position = position_jitter(width = .1,
                                        height = .1),
             size = 2,
             alpha = .4) +
  labs(x = "Puzzle order (1 through 8)",
       y = "% solved",
       title = "% solved over time",
       fill = "Puzzle") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1,8,1))

Fitting the mixed effects model:

# mixed effects model

lm.pct_solved_empty =
  glmer(solved ~
          1 +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = binomial(link = "logit"),
        data = df.analysis)

lm.pct_solved_order =
  glmer(solved ~
          1 +
          # fixed effect for puzzle order
          trialNum +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = binomial(link = "logit"),
        data = df.analysis)

# summary
lm.pct_solved_order %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## solved ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) +  
##     scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##    Data: df.analysis
## 
##      AIC      BIC   logLik deviance df.resid 
##   1840.4   1878.2   -913.2   1826.4     1633 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0237 -0.5676 -0.3125  0.6583  3.3794 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  subject_id (Intercept) 1.835    1.355   
## Number of obs: 1640, groups:  subject_id, 205
## 
## Fixed effects:
##                                   Estimate Std. Error z value
## (Intercept)                      -0.392330   0.168491  -2.328
## trialNum                         -0.098816   0.027951  -3.535
## scale(age)                       -0.218847   0.117493  -1.863
## scale(sokobanonline__solve_rate)  0.824983   0.071831  11.485
## scale(astar_iters)                0.004801   0.063904   0.075
## scale(participantEffort)          0.117835   0.115586   1.019
##                                  Pr(>|z|)    
## (Intercept)                      0.019886 *  
## trialNum                         0.000407 ***
## scale(age)                       0.062514 .  
## scale(sokobanonline__solve_rate)  < 2e-16 ***
## scale(astar_iters)               0.940118    
## scale(participantEffort)         0.307987    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g) s(___) scl(_)
## trialNum    -0.715                            
## scale(age)   0.015  0.017                     
## scl(skb___)  0.088 -0.241 -0.045              
## scl(str_tr)  0.090 -0.133 -0.005  0.144       
## scl(prtcpE) -0.007 -0.006 -0.118  0.024  0.004
anova(lm.pct_solved_empty, lm.pct_solved_order)
## Data: df.analysis
## Models:
## lm.pct_solved_empty: solved ~ 1 + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
## lm.pct_solved_order: solved ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##                     npar    AIC    BIC  logLik deviance  Chisq Df
## lm.pct_solved_empty    6 1851.1 1883.5 -919.54   1839.1          
## lm.pct_solved_order    7 1840.4 1878.2 -913.20   1826.4 12.691  1
##                     Pr(>Chisq)    
## lm.pct_solved_empty               
## lm.pct_solved_order  0.0003674 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate estimates
lm.pct_solved_order %>%
  ggpredict(bias_correction = T)
## $trialNum
## # Predicted probabilities of solved
## 
## trialNum | Predicted |     95% CI
## ---------------------------------
##        1 |      0.47 | 0.48, 0.51
##        2 |      0.46 | 0.48, 0.51
##        3 |      0.46 | 0.48, 0.50
##        4 |      0.45 | 0.48, 0.50
##        5 |      0.43 | 0.47, 0.50
##        6 |      0.42 | 0.47, 0.50
##        7 |      0.41 | 0.46, 0.49
##        8 |      0.39 | 0.45, 0.49
## 
## Adjusted for:
## *                       age =   36.74
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *         participantEffort =    4.80
## *                subject_id = 0 (population-level)
## 
## $age
## # Predicted probabilities of solved
## 
## age | Predicted |     95% CI
## ----------------------------
##  15 |      0.48 | 0.47, 0.52
##  20 |      0.47 | 0.47, 0.51
##  30 |      0.45 | 0.48, 0.50
##  40 |      0.43 | 0.47, 0.50
##  45 |      0.42 | 0.46, 0.50
##  50 |      0.41 | 0.45, 0.50
##  60 |      0.38 | 0.42, 0.50
##  75 |      0.34 | 0.35, 0.50
## 
## Adjusted for:
## *                  trialNum =    4.50
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *         participantEffort =    4.80
## *                subject_id = 0 (population-level)
## 
## $sokobanonline__solve_rate
## # Predicted probabilities of solved
## 
## sokobanonline__solve_rate | Predicted |     95% CI
## --------------------------------------------------
##                      0.40 |      0.03 | 0.04, 0.14
##                      0.60 |      0.21 | 0.31, 0.42
##                      0.80 |      0.50 | 0.49, 0.51
##                      1.00 |      0.76 | 0.54, 0.68
## 
## Adjusted for:
## *          trialNum =    4.50
## *               age =   36.74
## *       astar_iters = 6796.94
## * participantEffort =    4.80
## *        subject_id = 0 (population-level)
## 
## $astar_iters
## # Predicted probabilities of solved
## 
## astar_iters | Predicted |     95% CI
## ------------------------------------
##         650 |      0.44 | 0.47, 0.50
##        5500 |      0.44 | 0.47, 0.50
##       10300 |      0.44 | 0.47, 0.50
##       15150 |      0.44 | 0.47, 0.50
##       20000 |      0.44 | 0.47, 0.50
##       24800 |      0.44 | 0.47, 0.50
##       29650 |      0.44 | 0.46, 0.51
##       39300 |      0.44 | 0.45, 0.51
## 
## Adjusted for:
## *                  trialNum =  4.50
## *                       age = 36.74
## * sokobanonline__solve_rate =  0.73
## *         participantEffort =  4.80
## *                subject_id = 0 (population-level)
## 
## $participantEffort
## # Predicted probabilities of solved
## 
## participantEffort | Predicted |     95% CI
## ------------------------------------------
##                 1 |      0.30 | 0.17, 0.56
##                 3 |      0.38 | 0.37, 0.52
##                 4 |      0.42 | 0.45, 0.50
##                 5 |      0.44 | 0.48, 0.50
## 
## Adjusted for:
## *                  trialNum =    4.50
## *                       age =   36.74
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *                subject_id = 0 (population-level)
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "."

4.1.1.2 Boxes solved

plot2 =
  df.analysis %>%
  mutate(trialNum = case_when(trialNum > 4 ~ "Last half of puzzles",
                              trialNum <= 4 ~ "First half of puzzles"),
         puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                        solved,
                                        .fun = mean, .na_rm = TRUE)) %>%
  filter(!is.na(trialNum)) %>%
  ggplot(aes(x = puzzle_id_sorted,
             y = boxesSolved,
             color = trialNum)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "pointrange",
               position = position_dodge(.8)) +
  geom_point(position = position_jitter(width = .1,
                                        height = .1),
             size = 2,
             alpha = .1) +
  scale_y_continuous(limits = c(0, 3),
                     breaks = seq(0,3,1)) +
  labs(x = "Puzzle name",
       y = "# of boxes solved (0-3)",
       title = "# of boxes solved over time",
       color = "Stage") +
  coord_flip() +
  theme(legend.position = "top")

plot2

When broken down by individual puzzle:

df.analysis %>%
  group_by(trialNum, puzzle_id_sorted) %>%
  summarise(mean = mean(boxesSolved)) %>%
  ggplot(aes(x = trialNum,
             y = mean,
             color = puzzle_id_sorted)) +
  geom_point(position = position_jitter(width = .1,
                                        height = .1),
             size = 2,
             alpha = .4) +
  labs(x = "Puzzle order (1 through 8)",
       y = "Mean # of boxes solved",
       title = "Mean # of boxes solved over time",
       fill = "Puzzle") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1,8,1)) +
  scale_y_continuous(limits = c(0,3),
                     breaks = seq(0,3,1))

Fitting the mixed effects model:

# mixed effects model
lm.number_solved_empty =
  glmer(boxesSolved ~
          1 +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.analysis)

lm.number_solved_order =
  glmer(boxesSolved ~
          1 +
          # fixed effect for puzzle order (controlling for difficulty)
          trialNum +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.analysis)

# summary
lm.number_solved_order %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## boxesSolved ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) +  
##     scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##    Data: df.analysis
## 
##      AIC      BIC   logLik deviance df.resid 
##   4820.1   4857.9  -2403.1   4806.1     1633 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57643 -0.25869  0.02915  0.48384  1.06362 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  subject_id (Intercept) 0        0       
## Number of obs: 1640, groups:  subject_id, 205
## 
## Fixed effects:
##                                   Estimate Std. Error z value
## (Intercept)                       0.815209   0.037541  21.715
## trialNum                         -0.017653   0.007571  -2.332
## scale(age)                       -0.024591   0.017303  -1.421
## scale(sokobanonline__solve_rate)  0.087131   0.016979   5.132
## scale(astar_iters)               -0.037730   0.018226  -2.070
## scale(participantEffort)          0.010122   0.017409   0.581
##                                  Pr(>|z|)    
## (Intercept)                       < 2e-16 ***
## trialNum                           0.0197 *  
## scale(age)                         0.1553    
## scale(sokobanonline__solve_rate) 2.87e-07 ***
## scale(astar_iters)                 0.0384 *  
## scale(participantEffort)           0.5610    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g) s(___) scl(_)
## trialNum    -0.890                            
## scale(age)   0.011  0.000                     
## scl(skb___)  0.092 -0.144 -0.005              
## scl(str_tr)  0.127 -0.121 -0.004  0.106       
## scl(prtcpE) -0.001 -0.004 -0.107  0.005  0.005
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(lm.number_solved_empty, lm.number_solved_order)
## Data: df.analysis
## Models:
## lm.number_solved_empty: boxesSolved ~ 1 + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
## lm.number_solved_order: boxesSolved ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##                        npar    AIC    BIC  logLik deviance  Chisq Df
## lm.number_solved_empty    6 4823.6 4856.0 -2405.8   4811.6          
## lm.number_solved_order    7 4820.1 4857.9 -2403.1   4806.1 5.4415  1
##                        Pr(>Chisq)  
## lm.number_solved_empty             
## lm.number_solved_order    0.01966 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate estimates
lm.number_solved_order %>%
  ggpredict(bias_correction = T)
## $trialNum
## # Predicted counts of boxesSolved
## 
## trialNum | Predicted |     95% CI
## ---------------------------------
##        1 |      2.22 | 2.09, 2.36
##        2 |      2.18 | 2.08, 2.29
##        3 |      2.14 | 2.06, 2.23
##        4 |      2.11 | 2.03, 2.18
##        5 |      2.07 | 2.00, 2.14
##        6 |      2.03 | 1.95, 2.12
##        7 |      2.00 | 1.90, 2.10
##        8 |      1.96 | 1.84, 2.09
## 
## Adjusted for:
## *                       age =   36.74
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *         participantEffort =    4.80
## *                subject_id = 0 (population-level)
## 
## $age
## # Predicted counts of boxesSolved
## 
## age | Predicted |     95% CI
## ----------------------------
##  15 |      2.18 | 2.04, 2.33
##  20 |      2.16 | 2.04, 2.28
##  30 |      2.12 | 2.04, 2.20
##  40 |      2.07 | 2.00, 2.15
##  45 |      2.05 | 1.97, 2.14
##  50 |      2.03 | 1.93, 2.14
##  60 |      1.99 | 1.85, 2.14
##  75 |      1.93 | 1.73, 2.16
## 
## Adjusted for:
## *                  trialNum =    4.50
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *         participantEffort =    4.80
## *                subject_id = 0 (population-level)
## 
## $sokobanonline__solve_rate
## # Predicted counts of boxesSolved
## 
## sokobanonline__solve_rate | Predicted |     95% CI
## --------------------------------------------------
##                      0.40 |      1.44 | 1.24, 1.67
##                      0.60 |      1.80 | 1.68, 1.93
##                      0.80 |      2.26 | 2.16, 2.36
##                      1.00 |      2.83 | 2.51, 3.18
## 
## Adjusted for:
## *          trialNum =    4.50
## *               age =   36.74
## *       astar_iters = 6796.94
## * participantEffort =    4.80
## *        subject_id = 0 (population-level)
## 
## $astar_iters
## # Predicted counts of boxesSolved
## 
## astar_iters | Predicted |     95% CI
## ------------------------------------
##         650 |      2.14 | 2.05, 2.22
##        5500 |      2.10 | 2.03, 2.17
##       10300 |      2.06 | 1.98, 2.14
##       15150 |      2.02 | 1.93, 2.12
##       20000 |      1.98 | 1.87, 2.11
##       24800 |      1.95 | 1.80, 2.10
##       29650 |      1.91 | 1.74, 2.09
##       39300 |      1.84 | 1.62, 2.09
## 
## Adjusted for:
## *                  trialNum =  4.50
## *                       age = 36.74
## * sokobanonline__solve_rate =  0.73
## *         participantEffort =  4.80
## *                subject_id = 0 (population-level)
## 
## $participantEffort
## # Predicted counts of boxesSolved
## 
## participantEffort | Predicted |     95% CI
## ------------------------------------------
##                 1 |      1.94 | 1.50, 2.50
##                 3 |      2.01 | 1.78, 2.28
##                 4 |      2.05 | 1.93, 2.19
##                 5 |      2.10 | 2.02, 2.17
## 
## Adjusted for:
## *                  trialNum =    4.50
## *                       age =   36.74
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *                subject_id = 0 (population-level)
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "."

4.1.1.3 Time spent

plot3 =
  df.analysis %>%
  mutate(trialNum = case_when(trialNum > 4 ~ "Last half of puzzles",
                              trialNum <= 4 ~ "First half of puzzles"),
         puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                        solved,
                                        .fun = mean, .na_rm = TRUE)) %>%
  filter(!is.na(trialNum)) %>%
  ggplot(aes(x = puzzle_id_sorted,
             y = solveDuration,
             color = trialNum)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "pointrange",
               position = position_dodge(.8)) +
  geom_point(position = position_jitter(width = .1,
                                        height = .1),
             size = 2,
             alpha = .1) +
  scale_y_continuous(limits = c(0, 300),
                     breaks = seq(0,300,60)) +
  labs(x = "Puzzle name",
       y = "# of boxes solved (0-3)",
       title = "# of boxes solved over time",
       color = "Stage") +
  coord_flip() +
  theme(legend.position = "top")

plot3

When broken down by puzzle:

df.analysis %>%
  group_by(trialNum, puzzle_id_sorted) %>%
  summarise(mean = mean(solveDuration, na.rm = T)) %>%
  ggplot(aes(x = trialNum,
             y = mean,
             color = puzzle_id_sorted)) +
  geom_point(position = position_jitter(width = .1,
                                        height = .1),
             size = 2,
             alpha = .4) +
  labs(x = "Puzzle order (1 through 8)",
       y = "Mean solve time",
       title = "Mean solve time over time",
       fill = "Puzzle") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1,8,1)) +
  scale_y_continuous(limits = c(0,300),
                     breaks = seq(0,300,60))

Fitting the mixed effects model:

# mixed effects model
lm.duration_order_empty =
  glmer(solveDuration ~
          1 +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.analysis)

lm.duration_order =
  glmer(solveDuration ~
          1 +
          # fixed effect for puzzle order (controlling for difficulty)
          trialNum +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.analysis)

# summary
lm.duration_order %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## solveDuration ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) +  
##     scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##    Data: df.analysis
## 
##      AIC      BIC   logLik deviance df.resid 
##  19297.2  19327.9  -9641.6  19283.2      587 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -10.758  -3.800  -0.253   3.056  16.312 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  subject_id (Intercept) 0.1523   0.3903  
## Number of obs: 594, groups:  subject_id, 170
## 
## Fixed effects:
##                                   Estimate Std. Error z value
## (Intercept)                       4.964997   0.031484 157.700
## trialNum                          0.009484   0.001956   4.848
## scale(age)                        0.022585   0.031470   0.718
## scale(sokobanonline__solve_rate) -0.126157   0.004136 -30.505
## scale(astar_iters)                0.026036   0.004693   5.547
## scale(participantEffort)          0.011174   0.030190   0.370
##                                  Pr(>|z|)    
## (Intercept)                       < 2e-16 ***
## trialNum                         1.25e-06 ***
## scale(age)                          0.473    
## scale(sokobanonline__solve_rate)  < 2e-16 ***
## scale(astar_iters)               2.90e-08 ***
## scale(participantEffort)            0.711    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g) s(___) scl(_)
## trialNum    -0.264                            
## scale(age)   0.073 -0.007                     
## scl(skb___) -0.022 -0.091  0.002              
## scl(str_tr)  0.097 -0.356  0.005  0.098       
## scl(prtcpE) -0.017 -0.001 -0.151  0.004 -0.006
anova(lm.duration_order_empty, lm.duration_order)
## Data: df.analysis
## Models:
## lm.duration_order_empty: solveDuration ~ 1 + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
## lm.duration_order: solveDuration ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##                         npar   AIC   BIC  logLik deviance Chisq Df
## lm.duration_order_empty    6 19319 19345 -9653.3    19307         
## lm.duration_order          7 19297 19328 -9641.6    19283 23.49  1
##                         Pr(>Chisq)    
## lm.duration_order_empty               
## lm.duration_order        1.256e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate estimates
lm.duration_order %>%
  ggpredict(bias_correction = T)
## $trialNum
## # Predicted counts of solveDuration
## 
## trialNum | Predicted |         95% CI
## -------------------------------------
##        1 |    138.16 | 130.46, 147.28
##        2 |    139.48 | 131.79, 148.58
##        3 |    140.81 | 133.10, 149.93
##        4 |    142.15 | 134.40, 151.33
##        5 |    143.50 | 135.67, 152.78
##        6 |    144.87 | 136.92, 154.28
##        7 |    146.25 | 138.15, 155.84
##        8 |    147.64 | 139.36, 157.45
## 
## Adjusted for:
## *                       age =   35.78
## * sokobanonline__solve_rate =    0.76
## *               astar_iters = 6266.02
## *         participantEffort =    4.83
## *                subject_id = 0 (population-level)
## 
## $age
## # Predicted counts of solveDuration
## 
## age | Predicted |         95% CI
## --------------------------------
##  15 |    137.35 | 122.14, 155.47
##  20 |    138.62 | 125.89, 153.63
##  30 |    141.19 | 132.55, 151.37
##  40 |    143.81 | 135.51, 153.62
##  45 |    145.13 | 135.12, 156.91
##  50 |    146.47 | 134.01, 161.15
##  60 |    149.19 | 130.84, 171.24
##  75 |    153.36 | 125.36, 188.84
## 
## Adjusted for:
## *                  trialNum =    4.41
## * sokobanonline__solve_rate =    0.76
## *               astar_iters = 6266.02
## *         participantEffort =    4.83
## *                subject_id = 0 (population-level)
## 
## $sokobanonline__solve_rate
## # Predicted counts of solveDuration
## 
## sokobanonline__solve_rate | Predicted |         95% CI
## ------------------------------------------------------
##                      0.60 |    185.05 | 174.56, 197.45
##                      0.70 |    157.19 | 148.58, 167.40
##                      0.80 |    133.53 | 126.23, 142.18
##                      0.90 |    113.43 | 107.05, 120.98
## 
## Adjusted for:
## *          trialNum =    4.41
## *               age =   35.78
## *       astar_iters = 6266.02
## * participantEffort =    4.83
## *        subject_id = 0 (population-level)
## 
## $astar_iters
## # Predicted counts of solveDuration
## 
## astar_iters | Predicted |         95% CI
## ----------------------------------------
##         650 |    140.55 | 132.85, 149.67
##        5500 |    142.40 | 134.64, 151.60
##       10300 |    144.26 | 136.38, 153.59
##       15150 |    146.16 | 138.12, 155.69
##       20000 |    148.09 | 139.82, 157.87
##       24800 |    150.02 | 141.49, 160.11
##       29650 |    152.00 | 143.14, 162.46
##       39300 |    156.01 | 146.37, 167.37
## 
## Adjusted for:
## *                  trialNum =  4.41
## *                       age = 35.78
## * sokobanonline__solve_rate =  0.76
## *         participantEffort =  4.83
## *                subject_id = 0 (population-level)
## 
## $participantEffort
## # Predicted counts of solveDuration
## 
## participantEffort | Predicted |         95% CI
## ----------------------------------------------
##                 1 |    131.31 |  84.68, 204.96
##                 3 |    137.14 | 110.84, 170.80
##                 4 |    140.16 | 125.94, 157.00
##                 5 |    143.23 | 134.89, 153.09
## 
## Adjusted for:
## *                  trialNum =    4.41
## *                       age =   35.78
## * sokobanonline__solve_rate =    0.76
## *               astar_iters = 6266.02
## *                subject_id = 0 (population-level)
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "."

4.1.1.4 No of moves

plot4 =
  df.analysis %>%
  mutate(trialNum = case_when(trialNum > 4 ~ "Last half of puzzles",
                              trialNum <= 4 ~ "First half of puzzles"),
         puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                                        solved,
                                        .fun = mean, .na_rm = TRUE)) %>%
  filter(!is.na(trialNum)) %>%
  ggplot(aes(x = puzzle_id_sorted,
             y = attempt_nInputEvents,
             color = trialNum)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "pointrange",
               position = position_dodge(.8)) +
  geom_point(position = position_jitter(width = .1,
                                        height = .1),
             size = 2,
             alpha = .1) +
  labs(x = "Puzzle name",
       y = "# of moves",
       title = "# of moves over time",
       color = "Stage") +
  coord_flip() +
  theme(legend.position = "top")

plot4

combined_plot = 
  plot1 + plot2 + plot3 + plot4

combined_plot

When broken down by puzzle:

df.analysis %>%
  group_by(trialNum, puzzle_id_sorted) %>%
  summarise(mean = mean(attempt_nInputEvents, na.rm = T)) %>%
  ggplot(aes(x = trialNum,
             y = mean,
             color = puzzle_id_sorted)) +
  geom_point(position = position_jitter(width = .1,
                                        height = .1),
             size = 2,
             alpha = .4) +
  labs(x = "Puzzle order (1 through 8)",
       y = "Mean # of moves",
       title = "Mean # of moves over time",
       fill = "Puzzle") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1,8,1)) +
  scale_y_continuous(limits = c(0, 600),
                     breaks = seq(0,600,100))

Fitting the mixed effects model:

# mixed effects model
lm.moves_order_empty =
  glmer(attempt_nInputEvents ~
          1 +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.analysis)

lm.moves_order =
  glmer(attempt_nInputEvents ~
          1 +
          # fixed effect for puzzle order (controlling for difficulty)
          trialNum +
          # controlling for other variables
          scale(age) + # age
          scale(sokobanonline__solve_rate) + # difficulty
          scale(astar_iters) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.analysis)

# summary
lm.moves_order %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## attempt_nInputEvents ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) +  
##     scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##    Data: df.analysis
## 
##      AIC      BIC   logLik deviance df.resid 
##  61033.1  61070.9 -30509.6  61019.1     1633 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -14.4090  -3.7989  -0.2077   3.2561  24.7126 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  subject_id (Intercept) 0.2083   0.4563  
## Number of obs: 1640, groups:  subject_id, 205
## 
## Fixed effects:
##                                    Estimate Std. Error z value
## (Intercept)                       5.1546860  0.0321496 160.334
## trialNum                         -0.0206785  0.0008303 -24.904
## scale(age)                       -0.0645355  0.0321408  -2.008
## scale(sokobanonline__solve_rate) -0.1374895  0.0019654 -69.956
## scale(astar_iters)                0.0653442  0.0017487  37.366
## scale(participantEffort)         -0.0589342  0.0321208  -1.835
##                                  Pr(>|z|)    
## (Intercept)                        <2e-16 ***
## trialNum                           <2e-16 ***
## scale(age)                         0.0447 *  
## scale(sokobanonline__solve_rate)   <2e-16 ***
## scale(astar_iters)                 <2e-16 ***
## scale(participantEffort)           0.0665 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g) s(___) scl(_)
## trialNum    -0.114                            
## scale(age)   0.000  0.000                     
## scl(skb___)  0.029 -0.183 -0.001              
## scl(str_tr)  0.010 -0.121  0.000  0.087       
## scl(prtcpE)  0.000  0.000 -0.108  0.000  0.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
anova(lm.moves_order_empty, lm.moves_order)
## Data: df.analysis
## Models:
## lm.moves_order_empty: attempt_nInputEvents ~ 1 + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
## lm.moves_order: attempt_nInputEvents ~ 1 + trialNum + scale(age) + scale(sokobanonline__solve_rate) + scale(astar_iters) + scale(participantEffort) + (1 | subject_id)
##                      npar   AIC   BIC logLik deviance Chisq Df
## lm.moves_order_empty    6 61651 61684 -30820    61639         
## lm.moves_order          7 61033 61071 -30510    61019 620.3  1
##                      Pr(>Chisq)    
## lm.moves_order_empty               
## lm.moves_order        < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate estimates
lm.moves_order %>%
  ggpredict(bias_correction = T)
## $trialNum
## # Predicted counts of attempt_nInputEvents
## 
## trialNum | Predicted |         95% CI
## -------------------------------------
##        1 |    170.17 | 160.25, 181.71
##        2 |    166.69 | 156.99, 177.97
##        3 |    163.28 | 153.79, 174.32
##        4 |    159.93 | 150.65, 170.74
##        5 |    156.66 | 147.56, 167.25
##        6 |    153.45 | 144.54, 163.83
##        7 |    150.31 | 141.57, 160.50
##        8 |    147.24 | 138.65, 157.23
## 
## Adjusted for:
## *                       age =   36.74
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *         participantEffort =    4.80
## *                subject_id = 0 (population-level)
## 
## $age
## # Predicted counts of attempt_nInputEvents
## 
## age | Predicted |         95% CI
## --------------------------------
##  15 |    177.43 | 156.58, 202.18
##  20 |    172.83 | 155.85, 192.73
##  30 |    163.99 | 153.10, 176.63
##  40 |    155.60 | 146.24, 166.48
##  45 |    151.57 | 140.93, 163.93
##  50 |    147.64 | 134.98, 162.39
##  60 |    140.09 | 122.78, 160.73
##  75 |    129.48 | 105.68, 159.52
## 
## Adjusted for:
## *                  trialNum =    4.50
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *         participantEffort =    4.80
## *                subject_id = 0 (population-level)
## 
## $sokobanonline__solve_rate
## # Predicted counts of attempt_nInputEvents
## 
## sokobanonline__solve_rate | Predicted |         95% CI
## ------------------------------------------------------
##                      0.40 |    284.99 | 267.90, 304.86
##                      0.60 |    199.71 | 188.06, 213.27
##                      0.80 |    139.95 | 131.81, 149.43
##                      1.00 |     98.07 |  92.24, 104.86
## 
## Adjusted for:
## *          trialNum =    4.50
## *               age =   36.74
## *       astar_iters = 6796.94
## * participantEffort =    4.80
## *        subject_id = 0 (population-level)
## 
## $astar_iters
## # Predicted counts of attempt_nInputEvents
## 
## astar_iters | Predicted |         95% CI
## ----------------------------------------
##         650 |    151.83 | 143.01, 162.10
##        5500 |    156.90 | 147.79, 167.51
##       10300 |    162.09 | 152.68, 173.05
##       15150 |    167.51 | 157.77, 178.84
##       20000 |    173.11 | 163.03, 184.83
##       24800 |    178.83 | 168.40, 190.97
##       29650 |    184.81 | 173.99, 197.39
##       39300 |    197.30 | 185.65, 210.84
## 
## Adjusted for:
## *                  trialNum =  4.50
## *                       age = 36.74
## * sokobanonline__solve_rate =  0.73
## *         participantEffort =  4.80
## *                subject_id = 0 (population-level)
## 
## $participantEffort
## # Predicted counts of attempt_nInputEvents
## 
## participantEffort | Predicted |         95% CI
## ----------------------------------------------
##                 1 |    244.79 | 153.44, 392.73
##                 3 |    194.66 | 155.15, 245.59
##                 4 |    173.58 | 154.89, 195.62
##                 5 |    154.79 | 145.16, 165.98
## 
## Adjusted for:
## *                  trialNum =    4.50
## *                       age =   36.74
## * sokobanonline__solve_rate =    0.73
## *               astar_iters = 6796.94
## *                subject_id = 0 (population-level)
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "."

4.2 Success on preceding puzzle

Does success on the previous puzzle influence specific puzzle choice?

df.analysis =
  df.analysis %>%
  mutate(prev_puzzle_solved = lag(solved))

df.analysis %>%
  filter(!is.na(prev_puzzle_solved)) %>%
  mutate(puzzle_id_ordered = fct_reorder(puzzle_id_sorted, sokobanonline__solve_rate),
         prev_puzzle_solved = case_when(prev_puzzle_solved == 1 ~ "Yes",
                                        prev_puzzle_solved == 0 ~ "No")) %>%
  group_by(prev_puzzle_solved, stimuli_set, puzzle_id_sorted) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(prev_puzzle_solved, stimuli_set) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(x = factor(prev_puzzle_solved),
             y = (prop*100),
             fill = puzzle_id_sorted)) +
  geom_col(position = "dodge") +
  facet_wrap(~ stimuli_set + puzzle_id_sorted, ncol = 8) +
  labs(x = "Previous puzzle solved?",
       y = "% of times puzzle was chosen",
       fill = "Puzzle",
       title = "Player Puzzle Selection by Previous Puzzle Success") +
  theme(legend.position = "none")

df.analysis %>%
  mutate(puzzle_id_sorted = fct_reorder(puzzle_id_sorted,
                               solved,
                               .fun = mean, .na_rm = TRUE),
         prev_puzzle_solved = case_when(prev_puzzle_solved == 1 ~ "Yes",
                                        prev_puzzle_solved == 0 ~ "No")) %>%
  filter(!is.na(prev_puzzle_solved)) %>%
  ggplot(aes(x = puzzle_id_sorted,
             y = (solved * 100),
             fill = factor(prev_puzzle_solved))) +
  stat_summary(fun = mean,
               geom = "bar",
               width = .6,
               alpha = .5,
               position = position_dodge(.8)) +
  stat_summary(fun.data = mean_cl_boot,
               geom = "errorbar",
               position = position_dodge(.8)) +
  scale_fill_manual(values = c("Yes" = "darkgreen",
                               "No" = "darkred")) +
  labs(x = "Puzzle name",
       y = "% solved",
       fill = "Last puzzle solved?",
       title = "% solved by previous puzzle success",
       subtitle = "N = 205 subjects, 24 puzzles total (n = 67-70 each)") +
  coord_flip() +
  theme(legend.position = "right")

Fitting the mixed effects model:

# mixed effects model
df.subset = df.analysis %>%
  filter(!is.na(prev_puzzle_solved),
         !is.na(sokobanonline__solve_rate),
         !is.na(participantEffort))

lm.reduced =
  glmer(solved ~
          1 +
          scale(sokobanonline__solve_rate) +
          scale(participantEffort) +
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.subset)

lm.previous_performance =
  glmer(solved ~
          1 +
          # fixed effect for previous puzzle solved
          prev_puzzle_solved +
          # controlling for other variables
          scale(sokobanonline__solve_rate) + # difficulty
          scale(participantEffort) + # subjective effort
          # random intercept for player
          (1 | subject_id),
        family = poisson(link = "log"),
        data = df.subset)

anova(lm.reduced, lm.previous_performance, test = "Chisq")
## Data: df.subset
## Models:
## lm.reduced: solved ~ 1 + scale(sokobanonline__solve_rate) + scale(participantEffort) + (1 | subject_id)
## lm.previous_performance: solved ~ 1 + prev_puzzle_solved + scale(sokobanonline__solve_rate) + scale(participantEffort) + (1 | subject_id)
##                         npar    AIC    BIC  logLik deviance Chisq Df
## lm.reduced                 4 2294.8 2316.4 -1143.4   2286.8         
## lm.previous_performance    5 2282.4 2309.5 -1136.2   2272.4 14.39  1
##                         Pr(>Chisq)    
## lm.reduced                            
## lm.previous_performance  0.0001486 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summary
lm.previous_performance %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## solved ~ 1 + prev_puzzle_solved + scale(sokobanonline__solve_rate) +  
##     scale(participantEffort) + (1 | subject_id)
##    Data: df.subset
## 
##      AIC      BIC   logLik deviance df.resid 
##   2282.5   2309.5  -1136.2   2272.5     1634 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8948 -0.5286 -0.4356  0.6578  2.2635 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  subject_id (Intercept) 0.1656   0.407   
## Number of obs: 1639, groups:  subject_id, 205
## 
## Fixed effects:
##                                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)                      -1.30400    0.06778 -19.238  < 2e-16
## prev_puzzle_solved                0.33470    0.08726   3.836 0.000125
## scale(sokobanonline__solve_rate)  0.35504    0.03951   8.985  < 2e-16
## scale(participantEffort)          0.04056    0.05201   0.780 0.435429
##                                     
## (Intercept)                      ***
## prev_puzzle_solved               ***
## scale(sokobanonline__solve_rate) ***
## scale(participantEffort)            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) prv_p_ s(___)
## prv_pzzl_sl -0.502              
## scl(skb___) -0.226  0.048       
## scl(prtcpE) -0.014 -0.022 -0.007
# Calculate estimates
lm.previous_performance %>%
  ggpredict(bias_correction = T)
## $prev_puzzle_solved
## # Predicted counts of solved
## 
## prev_puzzle_solved | Predicted |     95% CI
## -------------------------------------------
##                  0 |      0.45 | 0.66, 0.85
##                  1 |      0.63 | 0.91, 1.21
## 
## Adjusted for:
## * sokobanonline__solve_rate = 0.73
## *         participantEffort = 4.80
## *                subject_id = 0 (population-level)
## 
## $sokobanonline__solve_rate
## # Predicted counts of solved
## 
## sokobanonline__solve_rate | Predicted |     95% CI
## --------------------------------------------------
##                      0.40 |      0.11 | 0.13, 0.27
##                      0.60 |      0.28 | 0.38, 0.56
##                      0.80 |      0.70 | 1.04, 1.30
##                      1.00 |      1.75 | 2.24, 3.80
## 
## Adjusted for:
## * prev_puzzle_solved = 0.36
## *  participantEffort = 4.80
## *         subject_id = 0 (population-level)
## 
## $participantEffort
## # Predicted counts of solved
## 
## participantEffort | Predicted |     95% CI
## ------------------------------------------
##                 1 |      0.38 | 0.29, 1.35
##                 3 |      0.44 | 0.50, 1.07
##                 4 |      0.48 | 0.65, 0.96
##                 5 |      0.52 | 0.77, 0.96
## 
## Adjusted for:
## *        prev_puzzle_solved = 0.36
## * sokobanonline__solve_rate = 0.73
## *                subject_id = 0 (population-level)
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "."

4.3 Growth curve model

% solved:

# null model: no learning over time
lm.solve_growth_null =
  glmer(solved ~
          1 +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
        data = df.analysis,
        family = binomial(link = "logit"),
        control = glmerControl(optimizer = "bobyqa")) 

lm.solve_growth =
  glmer(solved ~
          1 +
          # fixed effect
          trialNum +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
        data = df.analysis,
        family = binomial(link = "logit"),
        control = glmerControl(optimizer = "bobyqa"))

lm.solve_growth %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## solved ~ 1 + trialNum + scale(age) + scale(participantEffort) +  
##     (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##    Data: df.analysis
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1544.0   1587.2   -764.0   1528.0     1632 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3715 -0.4343 -0.1516  0.3906  8.6867 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev. Corr
##  subject_id       (Intercept) 1.77503  1.3323       
##                   trialNum    0.01016  0.1008   1.00
##  puzzle_id_sorted (Intercept) 4.25909  2.0638       
## Number of obs: 1640, groups:  subject_id, 205; puzzle_id_sorted, 24
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)  
## (Intercept)              -1.183170   0.476426  -2.483   0.0130 *
## trialNum                 -0.008627   0.041586  -0.207   0.8357  
## scale(age)               -0.304039   0.143353  -2.121   0.0339 *
## scale(participantEffort)  0.136793   0.141100   0.969   0.3323  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g)
## trialNum    -0.340              
## scale(age)   0.015 -0.012       
## scl(prtcpE) -0.006  0.012 -0.119
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(lm.solve_growth_null, lm.solve_growth, test = "Chisq")
## Data: df.analysis
## Models:
## lm.solve_growth_null: solved ~ 1 + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
## lm.solve_growth: solved ~ 1 + trialNum + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##                      npar  AIC    BIC  logLik deviance Chisq Df
## lm.solve_growth_null    7 1542 1579.9 -764.02     1528         
## lm.solve_growth         8 1544 1587.2 -764.00     1528 0.043  1
##                      Pr(>Chisq)
## lm.solve_growth_null           
## lm.solve_growth          0.8357
df.analysis$predicted_solve =
  predict(lm.solve_growth, type = "response")

df.analysis_subject =
  df.analysis %>%
  group_by(subject_id, trialNum) %>%
  summarise(mean_pred = mean(predicted_solve, na.rm = TRUE),
            .groups = "drop")

plot1 =
  df.analysis_subject %>%
  ggplot(aes(x = trialNum,
             y = (mean_pred*100))) +
  geom_line(aes(group = subject_id),
            alpha = 0.1, color = "gray") +  # individual trajectories
  stat_summary(fun = mean,
               geom = "line",
               size = 1,
               color = "#1f77b4") +  # average curve
  stat_summary(fun.data = mean_cl_boot,
               geom = "ribbon",
               aes(group = 1),
               fill = "#1f77b4",
               alpha = 0.2) +
  labs(x = "Puzzle #",
       y = "% solved",
       title = "Predicted performance curves across puzzles: % solved",
       subtitle = "Grey lines are predicted lines for individual subjects,\nblue line and ribbon equal predicted mean + 95% CI") +
  scale_x_continuous(limits = c(1, 8),
                     breaks = seq(1,8,1)) +
  theme_minimal()

plot1

boxes solved:

# null model: no learning over time
lm.boxes_growth_null =
  glmer(boxesSolved ~
          1 +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
  data = df.analysis,
  family = poisson(link = "log"),
  control = glmerControl(optimizer = "bobyqa"))# full model: 

lm.boxes_growth =
  glmer(boxesSolved ~
          1 +
          trialNum +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
        data = df.analysis,
        family = poisson(link = "log"),
        control = glmerControl(optimizer = "bobyqa"))

lm.boxes_growth %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## boxesSolved ~ 1 + trialNum + scale(age) + scale(participantEffort) +  
##     (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##    Data: df.analysis
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   4759.0   4802.2  -2371.5   4743.0     1632 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.58261 -0.24194 -0.00553  0.33846  1.62867 
## 
## Random effects:
##  Groups           Name        Variance  Std.Dev.  Corr
##  subject_id       (Intercept) 0.000e+00 0.000e+00     
##                   trialNum    4.501e-16 2.122e-08  NaN
##  puzzle_id_sorted (Intercept) 5.070e-02 2.252e-01     
## Number of obs: 1640, groups:  subject_id, 205; puzzle_id_sorted, 24
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               0.727597   0.063852  11.395   <2e-16 ***
## trialNum                 -0.002063   0.009041  -0.228    0.820    
## scale(age)               -0.025980   0.017257  -1.505    0.132    
## scale(participantEffort)  0.008482   0.017548   0.483    0.629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g)
## trialNum    -0.638              
## scale(age)   0.008 -0.002       
## scl(prtcpE)  0.001 -0.005 -0.112
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(lm.boxes_growth_null, lm.boxes_growth, test = "Chisq")
## Data: df.analysis
## Models:
## lm.boxes_growth_null: boxesSolved ~ 1 + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
## lm.boxes_growth: boxesSolved ~ 1 + trialNum + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##                      npar  AIC    BIC  logLik deviance  Chisq Df
## lm.boxes_growth_null    7 4757 4794.8 -2371.5     4743          
## lm.boxes_growth         8 4759 4802.2 -2371.5     4743 0.0519  1
##                      Pr(>Chisq)
## lm.boxes_growth_null           
## lm.boxes_growth          0.8199
df.analysis$predicted_boxes =
  predict(lm.boxes_growth, type = "response")

df.analysis_subject =
  df.analysis %>%
  group_by(subject_id, trialNum) %>%
  summarise(mean_pred = mean(predicted_boxes, na.rm = TRUE),
            .groups = "drop")

plot2 =
  df.analysis_subject %>%
  ggplot(aes(x = trialNum,
             y = mean_pred)) +
  geom_line(aes(group = subject_id),
            alpha = 0.1, color = "gray") +  # individual trajectories
  stat_summary(fun = mean,
               geom = "line",
               size = 1,
               color = "#1f77b4") +  # average curve
  stat_summary(fun = "errorbar",
               geom = "line",
               size = 1,
               color = "#1f77b4") +  # average curve
  stat_summary(fun.data = mean_cl_boot,
               geom = "ribbon",
               aes(group = 1),
               fill = "#1f77b4",
               alpha = 0.2) +
  labs(x = "Puzzle #",
       y = "# of boxes solved",
       title = "Predicted performance curves across puzzles: # of boxes solved",
       subtitle = "Grey lines are predicted lines for individual subjects,\nblue line and ribbon equal predicted mean + 95% CI") +
  scale_x_continuous(limits = c(1, 8),
                     breaks = seq(1,8,1)) +
  scale_y_continuous(limits = c(0, 3),
                     breaks = seq(0,3, 1)) +
  theme_minimal()

plot2

time to solution:

# null model: no learning over time
lm.time_growth_null =
  glmer(solveDuration ~
          1 +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
  data = df.analysis,
  family = poisson(link = "log"),
  control = glmerControl(optimizer = "bobyqa"))# full model: 

lm.time_growth =
  glmer(solveDuration ~
          1 +
          trialNum +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
        data = df.analysis,
        family = poisson(link = "log"),
        control = glmerControl(optimizer = "bobyqa"))

lm.time_growth %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## solveDuration ~ 1 + trialNum + scale(age) + scale(participantEffort) +  
##     (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##    Data: df.analysis
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  13095.4  13130.5  -6539.7  13079.4      586 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.7278 -2.2725 -0.0664  1.6040 12.9123 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev. Corr 
##  subject_id       (Intercept) 0.83535  0.9140        
##                   trialNum    0.03517  0.1875   -0.87
##  puzzle_id_sorted (Intercept) 0.11106  0.3333        
## Number of obs: 594, groups:  subject_id, 170; puzzle_id_sorted, 23
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               5.246034   0.102840  51.011   <2e-16 ***
## trialNum                 -0.040481   0.015905  -2.545   0.0109 *  
## scale(age)                0.005152   0.038299   0.135   0.8930    
## scale(participantEffort)  0.021529   0.035807   0.601   0.5477    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g)
## trialNum    -0.643              
## scale(age)   0.023  0.000       
## scl(prtcpE)  0.007 -0.005 -0.152
anova(lm.time_growth_null, lm.time_growth, test = "Chisq")
## Data: df.analysis
## Models:
## lm.time_growth_null: solveDuration ~ 1 + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
## lm.time_growth: solveDuration ~ 1 + trialNum + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##                     npar   AIC   BIC  logLik deviance Chisq Df
## lm.time_growth_null    7 13100 13130 -6542.9    13086         
## lm.time_growth         8 13095 13130 -6539.7    13079 6.424  1
##                     Pr(>Chisq)  
## lm.time_growth_null             
## lm.time_growth         0.01126 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df.analysis =
  df.analysis %>%
  mutate(predicted_time = NA_real_)

df.analysis$predicted_time[df.analysis$solved == 1] =
  predict(lm.time_growth, type = "response")

df.analysis_subject =
  df.analysis %>%
  group_by(subject_id, trialNum) %>%
  summarise(mean_pred = mean(predicted_time, na.rm = TRUE),
            .groups = "drop")

plot3 =
  df.analysis_subject %>%
  ggplot(aes(x = trialNum,
             y = mean_pred)) +
  geom_line(aes(group = subject_id),
            alpha = 0.1, color = "gray") +  # individual trajectories
  stat_summary(fun = mean,
               geom = "line",
               size = 1,
               color = "#1f77b4") +  # average curve
  stat_summary(fun = "errorbar",
               geom = "line",
               size = 1,
               color = "#1f77b4") +  # average curve
  stat_summary(fun.data = mean_cl_boot,
               geom = "ribbon",
               aes(group = 1),
               fill = "#1f77b4",
               alpha = 0.2) +
  labs(x = "Puzzle #",
       y = "Time until solution",
       title = "Predicted performance curves across puzzles: time until solution",
       subtitle = "Grey lines are predicted lines for individual subjects,\nblue line and ribbon equal predicted mean + 95% CI") +
  scale_x_continuous(limits = c(1, 8),
                     breaks = seq(1,8,1)) +
  scale_y_continuous(limits = c(0, 300),
                     breaks = seq(0, 300, 60)) +
  theme_minimal()

plot3

number of moves:

# null model: no learning over time
lm.moves_growth_null =
  glmer(attempt_nInputEvents ~
          1 +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
  data = df.analysis,
  family = poisson(link = "log"),
  control = glmerControl(optimizer = "bobyqa"))# full model: 

lm.moves_growth =
  glmer(attempt_nInputEvents ~
          1 +
          trialNum +
          # control vars
          scale(age) + # age
          scale(participantEffort) + # subjective effort
          (1 | puzzle_id_sorted) + # random slope for puzzle
          (1 + trialNum | subject_id), # random slopes and intercepts per subject
        data = df.analysis,
        family = poisson(link = "log"),
        control = glmerControl(optimizer = "bobyqa"))

lm.moves_growth %>%
  summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## attempt_nInputEvents ~ 1 + trialNum + scale(age) + scale(participantEffort) +  
##     (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##    Data: df.analysis
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  44561.1  44604.4 -22272.6  44545.1     1632 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -14.8539  -2.8979  -0.1722   2.5324  21.2875 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev. Corr 
##  subject_id       (Intercept) 0.26276  0.51260       
##                   trialNum    0.00760  0.08718  -0.47
##  puzzle_id_sorted (Intercept) 0.05728  0.23934       
## Number of obs: 1640, groups:  subject_id, 205; puzzle_id_sorted, 24
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               5.263712   0.060814  86.554   <2e-16 ***
## trialNum                 -0.053089   0.006198  -8.565   <2e-16 ***
## scale(age)               -0.075497   0.032157  -2.348   0.0189 *  
## scale(participantEffort) -0.062008   0.031997  -1.938   0.0526 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trilNm scl(g)
## trialNum    -0.287              
## scale(age)  -0.001  0.001       
## scl(prtcpE)  0.000  0.000 -0.108
anova(lm.moves_growth_null, lm.moves_growth, test = "Chisq")
## Data: df.analysis
## Models:
## lm.moves_growth_null: attempt_nInputEvents ~ 1 + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
## lm.moves_growth: attempt_nInputEvents ~ 1 + trialNum + scale(age) + scale(participantEffort) + (1 | puzzle_id_sorted) + (1 + trialNum | subject_id)
##                      npar   AIC   BIC logLik deviance  Chisq Df
## lm.moves_growth_null    7 44622 44660 -22304    44608          
## lm.moves_growth         8 44561 44604 -22273    44545 62.922  1
##                      Pr(>Chisq)    
## lm.moves_growth_null               
## lm.moves_growth        2.15e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df.analysis$predicted_moves =
  predict(lm.moves_growth, type = "response")

df.analysis_subject =
  df.analysis %>%
  group_by(subject_id, trialNum) %>%
  summarise(mean_pred = mean(predicted_moves, na.rm = TRUE),
            .groups = "drop")

plot4 =
  df.analysis_subject %>%
  ggplot(aes(x = trialNum,
             y = mean_pred)) +
  geom_line(aes(group = subject_id),
            alpha = 0.1, color = "gray") +  # individual trajectories
  stat_summary(fun = mean,
               geom = "line",
               size = 1,
               color = "#1f77b4") +  # average curve
  stat_summary(fun = "errorbar",
               geom = "line",
               size = 1,
               color = "#1f77b4") +  # average curve
  stat_summary(fun.data = mean_cl_boot,
               geom = "ribbon",
               aes(group = 1),
               fill = "#1f77b4",
               alpha = 0.2) +
  labs(x = "Puzzle #",
       y = "# of moves",
       title = "Predicted performance curves across puzzles: # of moves",
subtitle = "Grey lines are predicted lines for individual subjects,\nblue line and ribbon equal predicted mean + 95% CI") +
  scale_x_continuous(limits = c(1, 8),
                     breaks = seq(1,8,1)) +
  scale_y_continuous(limits = c(0, 600),
                     breaks = seq(0, 600, 100)) +
  theme_minimal()

plot4

combined_plot =
  plot1 + plot2 + plot3 + plot4

combined_plot

5